## HIT 2007 Reanalysis

library(RItools)
source("code/seqblock1.R")
source("code/seqblock2k.R")
source("code/mahal.R")

x <- read.table("data/horimatan07.txt", header = TRUE)
## trim to units of interest:
x <- x[x$D == x$Z,]     ## assigned == received
x <- x[!(is.na(x$Y)),]  ## Y observed
x$id <- 1:nrow(x)

d2.p <- xBalance(Z ~ block + age + college + past.turnout, data = x, report = c("chisquare.test"))$overall$p.value

## Sequentially block HIT data (optional -- loaded below):
d2.p.HIT <- array(NA, 100)
for(i in 1:100){
  print(i)
  ## given dat X, do an SB.      
  sb1 <- seqblock1(id.vars = "id", id.vals = x[1, "id"], covar.vars = c("block", "age", "college", "past.turnout"), covar.vals = x[1, c("block", "age", "college", "past.turnout")], tr.names = c("0","1"), assg.prob.kfac = 5, file.name = "~/Desktop/hitTest.RData", verbose = FALSE)
  
  for(n.idx in 2:nrow(x)){
sb2k <- seqblock2k(object = "~/Desktop/hitTest.RData", id.vals = x[n.idx, "id"],  covar.vals = x[n.idx, c("block", "age", "college", "past.turnout")], file.name = "~/Desktop/hitTest.RData", verbose = FALSE)
  }

  d2.p.HIT[i] <- xBalance(I(as.numeric(Tr)) ~ block + age + college + past.turnout, data = sb2k$orig, report = c("chisquare.test"))$overall$p.value
}

load("data/d2pHIT100.RData")

opar <- par()

## Figure 7 (left, color):
plot((1:100)/101, as.numeric(sort(d2.p.HIT)), xlim=c(0,1), ylim=c(0,1), col="blue", pch=18, ylab = bquote(d^2 ~ "p-values"), xlab="Sequential Reallocation", axes=FALSE, main = "Policy Information Experiment")
axis(1, at=c(1, 20, 40, 60, 80, 100)/101, labels=c(1, 20, 40, 60, 80, 100))
axis(2)
abline(h=d2.p, col="red")

## Figure 7 (left, black & white):
par(mar=c(5.1, 4.2, 4.1, 2.1))
plot((1:100)/101, as.numeric(sort(d2.p.HIT)), xlim=c(0,1), ylim=c(0,1), pch=18, ylab = bquote(d^2 ~ "p-values"), xlab="Sequential Reallocation", axes=FALSE, main = "Policy Information Experiment")
axis(1, at=c(1, 20, 40, 60, 80, 100)/101, labels=c(1, 20, 40, 60, 80, 100))
axis(2)
segments(x0 = 1/101, x1 = 100/101, y0 = d2.p, col = "grey", lwd = 2)

par(opar)
